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SUMMARY 


The propagation of a laser beam in an optically dense medium such 
as a fog, dust storm, or smoke, is a problem of growing importance, both 
for communication and detection purposes. Although such dense media 
lead to a significant attenuation of the primary beam, much of the 
scattered radiation may still be found close to the beam axis and will, 
thus, be available for detection by a suitable detector. In this report 
we examine the spreading of a laser beam using the small- angle scattering 
approximation to the equation of transfer. This approximation, which 
assumes that most photons travel essentially parallel to the beam axis, 
has also been used to study the propagation of fast charged particles 
through metal foils. It appears to be equally suited to the study of 
the propagation of beams of visible or near infrared light through media 
such as fog, dust or smoke, where the scattering phase function is 
highly anisotropic. 

This equation of transfer may be solved in closed form by the use 
of Fourier transform techniques. The resulting expressions are simplest 
for the radiance, or alternately the power received by a coaxial detector, 
rather than for the irradiance. In view of the assumptions involved in 
the small-angle approximation, it is the radiance which is of more 
interest anyway. The resulting expression involves a single integral 
from zero to infinity; the appendix outlines the procedure for its 
numerical evaluation. 


V 


In keeping with the approximate nature of our solution , and in order 
to fully exploit its mathematical simplicity, we have chosen simple analytic 
models for the forward peak of the scattering phase function, rather than 
using Mie theory. As well as the well-known Gaussian functional form, we 
also examine some exponential and binomial models. Our computational 
results indicate that , provided the parameters of the models are suitably 
selected, there is little to choose between the models, with the exception 
of the sea-water model, which we do not recommend. 

Despite the relative simplicity of the expressions that we obtain, a 
number of authors have resorted to further approximations, in order to 
extract even simpler results. The method of Dolin and Fante starts by 
separating the scattered and unscattered beams. In the case of a 
Gaussian phase function, this method leads to a single finite integral, 
which shows reasonable agreement with our results. The method of 
Arnush and Stotts is essentially a low-frequency approximation, which 
yields reasonable results well away from the beam axis but leads to 
unphysical results close to the axis. The method of Tam and Zardecki 
involves a series expansion of our integral leading to a series of 
multidimensional finite integrals ; it is applicable to both radiance 
and irradiance {which is its main advantage) , though only for the 
Gaussian phase function. 


VI 



1.0 INTRODUCTION 


If a relatively narrow beam propagates in a scattering mediiim, photons 
are constantly removed from the beam. However, if the scatterers are of 
a size equal to or greater than the radiation wavelength, such as in the 
case of smoke, dust or fog particles compared to visible wavelengths, then 
most scattering events will result in a comparatively small deflection of 
the photon. This may lead to a gradual spreading of the original beam, both 
in thickness and angle. 

In this report we examine the broadening of a laser beam , and the signal 
that may be detected, as functions of both experimental geometry and the 
properties of the scattering medium. We shall employ the small-angle approxi- 
mation to the equation of radiative transfer, which ignores photons which 
have suffered large deflections as they will be assumed lost. In order to 
obtain tractable answers, it will prove necessary to assume simple analytic 
forms for the scattering phase function and the initial beam profile. Never- 
theless, the analysis presented in this report will be as free as possible 
of unnecessary approximations. 


2.0 EQUATION OF TRANSFER IN SMALL-ANGLE APPROXIMATION 
Let I(z, r, fi) dy dn be the intensity of radiation (or the number of 
photons) in a volume element dV centered at the point z, r = (x, y) , and 
travelling within a cone of solid angle dn centered about the direction n. 
Then I satisfies the equation of radiative transfer, which we may write 
(Refs, land 2) I 


O I(z, r, n) = w^cr j P(n . n') I(z, r, n') dn' 


( 1 ) 


where o is the extinction coefficient (km ) 


2 


is the albedo of single scattering 
= cos ^ (n . n') is the scattering angle 
and P(4 j) is the scattering phase function 

Even allowing for cylindrical symmetry about the axis of propagation 
(the z axis ) , Eq. (1) is exceedingly hard to solve, even numerically. However, 
since the diameter of our detector will always be small compared to the total 
propagation distance, we may safely assume that all photons which are eventually 
detected will have spent their flight time travelling essentially parallel to 
the z-axis . We may thus set cos 0=1, where 0 is the angle the photon makes 
with the z axis . 

Note that this approximation ignores the contribution from all photons 
which vmdergo at least one large-angle scattering event. All such 
photons will clearly need to undergo at least a second large-angle scattering 
event, and maybe even a third, in order for them to reach the detector. As 
we are assuming that the phase function, P, is strongly forward-peaked, the 
probability of two or more large-angle scatterings is clearly very small, 
and thus the neglected contribution will be small. 

The main effect of this assumption is to replace the unit propagation 
vector by 

n - ( 1 ^ 2 ' * (2a) 

Although this new propagation vector is no longer correctly normalized, this 
should not cause any problems, as the number of photons for which in | 
is not true will clearly be small. 


« 1 
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The second effect is that we may use nj^ - as the argument of P 
in Eq. (1) , i. e. , 

P(n -n') P(n - n') . (2b) 

^ ~ ^ -_j_ 

The third effect is to replace the limits of this (two-dimensional) integral 
by + With these points in mind, we may now rewrite Eq. (1) as 



n 





I(z, 




n ) 


to 


o 


00 

f 

a 

. 

—oo 



n') I(z, r, n')dn' . 

U ~ -^1 


(3) 


Equation (3) is referred to as the equation of radiative transfer in the 
small-angle approximation. Its main advantage over Eq. (1) is in the simplification 
of the directional derivative. This equation has been used extensively in the 
theory of foil penetration by fast charged particles (Refs. 3, 4 and 5) . Though 
Wentzel (Ref. 3) was the first to use the small-angle approximation for charged 
particle transfer, perhaps the first person to employ this equation in the field 
of radiative transfer appears to be Dolin (Ref. 2) . 

One further result of the small-angle approximation is that all detected 
photons are assumed to have travelled the same distance. Thus their time of 
travel is constant, and a pulse will undergo no time-dispersion. 


3.0 FORMAL SOLUTION IN THE SMALL-ANGLE APPROXIMATION 
Equation (3) may be solved, at least formally, by the use of Fourier 
transform techniques. Introducing the definitions 


II nil I II I II 


I(z,n, O = ( 2 tt) 


-2 


T/ ^ + S*n.) , , 

I(z, r, rij^) e ~ ~ dr drij^ 


(4) 


-1 


fCO 


i ^-n. 

P (n^) e ~ ~J_ ^ 


and P(^) = (27T) 

we take the double Fourier transform of Eq. (3) to obtain (Ref. 2) 


(5) 


3 3 1^ /s /s 

^ - n . + crj I(z, ri/ = 27 t co^a P (^) l(z, - 


( 6 ) 


Equation (6) is easily solved, to yield (Refs. 2, 5) 


l(z, r)/0 = I (ri/ C + z n) e 
^ o ~ ~ ~ 


- OZ + ^ 


(7) 


where 


rz ^ 

^ = Q n) = 2tt 0) a P (|5 + z' ril) dz' 

o Jo ~ 


( 8 ) 


and 1^(13/ O is the Fourier transform of the initial intensity distribution 
(incident beam profile) at z = 0. 

To obtain the intensity distribution at any point in the medium, it is 
merely necessary to re-transform Eq. (7) (Refs. 6 and 7) 


-2 

► 

» 

■ f 

r, n ) = (2 tt) 

^ 'X-L « 





- i (n . r + 4 ■ n ) 

I(z, n# 5) e ~ ~ _L ^0 . (9) 
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The principal difficulty with this procedure is the evaluation, and 
subsequent exponentiation, of the function Evaluation of in the case of 
a real (Mie) phase function would appear to be prohibitive, and we shall employ 
instead a nimber of simpler analytic functions for P (described in later sections) . 
However, even with such simplifying assumptions, many authors have still founc. 
it necessary to employ further approximations in order to obtain tractable 
expressions for I (Eq. 7) . We shall examine two of these approximations later. 
Before we proceed to specific examples, however, we remind ourselves that 
it is irradiance (flux density) and received power, rather than radiance, which 
is of concern to us in this study. Using the relation between irradiance, N, 
and radiance, I, we may simplify Eq. (9) (Refs. 2 and 6) . 


N(z, r) 


l(z, r, 


27T 


n^) [n 


z] dn 


00 

‘ l(z, 

~oo 


r. 


n ) 



( 10 ) 


i(z, n, 0) 


-1 ri»r , 
e ~ ~ dn 


( 11 ) 


With the elimination of C in Eq. (11) , we may simplify Eqs. (7) and (8): 

^ -Gz + 

i(z, n, 0) = I (n, z n) e ° 

~ ~ o ~ ~ 


( 7 ’) 
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where 


= 27Twa f P(|z* Ti|) dz' 

o Jq 


Equation (11) can now be further simplified by an appeal to symmetry. Since 
P(\(;) clearly depends only on the scalar ]n |/ and not on the vector n , P will 
similarly be a scalar function, as will Similarly, if we assume that the 

incident beam profile is circularly symmetric, then 1^ (h / z Tl) will also be 
a scalar function of n. Thus I(z, r\, 0) will be a scalar function, and Eq. (11) 


becomes 


N(z, r) = 2 tt J (n r) I(z, n. 0) n dri . (11') 

Jq ° 

From Eq. (11') we see immediately that N is a scalar function of r, as we 

would expect from the above symmetry arguments. A more tractable expression 

for the case of a 5-function beam will be given in Eq. (31) . 

In most instances, of course, what we are most interested in (and what 

we physically measure) is the power received by some detector. This will 

involve the integration of Eq. (11') over the area of the detector, perhaps 

modulated by a response function. If we assume a coaxial, circular detector 

of radius R, with a flat response, then we have 

|-R 

P(z, R) = 2 tt N(z, r) r dr 

•'O 

= 47t^ e“^^ J, (n R) I (n, z n) e Or dn . ( 12 ) 

1 o 


0 
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This result should prove amenable to numerical integration, especially if 
a relatively simple expression for can be obtained. Sample results will 
be presented below. 

One useful result which can be obtained analytically, is the total power 
crossing a surface z = constant: 


P(z, oo) = 


N(z, r) dr 




• r 



. 


^ ^ -i n • r , 

I(z, ri, 0)e ~ ~dr|dr 


2 ^ 

= 4tt I(z, 0, 0) 


^ -az + 2 tt 00 a z p ( o) 

2 ^ o 

= 4tt I ( 0, 0 ) e 
o 


-( 1-00 ) az 

= F e ° (13) 

o 

^ -1 

where F is the incident total power, and we have used the fact that P(0) = (27T) 
o 

From Eq. (13) we see that the only energy removed from the beam is that lost 
by absorption — i.e., there is no backscatter. 

One further parameter which will often prove useful is the beam spread, 
which we may define as 


<r^> = 


N(z, r) r dr / 


N(z, r) r dr 


•'O 


, ( 1-00 ) az 

= 2Tr F e ° 

° -0 


N(z, r) r dr . 


(14) 
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4.0 EXACT SOLUTION FOR GAUSSIAN BEAMS 
At the entrance to a scattering medium, a laser beam profile can often 
be adequately represented by a Gaussian functional form, both for the radial 
distribution, and the angular divergence. Thus we have 

I (r, n ) = F 3^ exp(-3^ 9? - r^) . (15) 

o ~ ~J_ o -L 

This may be easily transformed, and, in particular, we have 


O O O O O O 

I (n, z n) = F ( 277 )" exp(-n / 4 y - z n/43). (i6) 

o o 

In general, the laser beam profile will be well collimated, so that 3 and Y 
will be large. The inclusion of Eq. (16) in Eq. (12) will in no way complicate 
the numerical integration, though in our examples later we will allow both to 
go to infinity, so as to reduce the number of parameters whose influence should 
be examined. In all practical calculations, however, realistic values of both 
parameters should be included. 

4.1 Gaussian Phase Functions 

A Gaussian functional form is also often employed to describe P(lJj) , since 
exact Mie theory is clearly somewhat impractical. Thus we choose to write 


P(ilj) 


2 a 


2,2 

-a ^ 


/ 2tt . 


(17) 
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The a is an adjustable parameter, which controls the shape of the forward peak, 

2 " 

and is related to the rms scattering angle \p defined as 



2tt 


P(i|;) ip dip 




(17a) 


It is easily shown that for the Gaussian case 




2 


-2 

a 


(17b) 


(Though a will usually be large, it will rarely, if ever, be as large as 
3 or y.) Taking the Fourier transform of Eq. (17), we find 


P(5) 


.00 

J (C i|J) P(4^) ip dip 


2 

e"^ /4a ^ 2 ^ 


(5M 


(18) 


and hence = o) CJ n ^ a /tT erf (zri/2a) 

o o 

where erf is the well-known error function. 
Substituting Eqs. (16) and (19) into (12) we find 


p(z, R) = F R J, (n R) exp foo O CX /iT erf (zn/2a) 

o Jo 1 o 

- 0z - ri^/4 - z^ ri^/4 3^] dri . 


(19) 


( 20 ) 
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Equation (20) is an exact equation, which is solved nxunerically for various 
values of B and y* Here, we make a simplification in Eq. (20) and consider 
the limiting case of 3, y This physically implies that the beam is 

collimated and has zero width in space. Then by making the variable changes 

T = a z 


T = Oi T 

s o 


G 


R/z 



X = n R 

Equation (20) can be reduced to 


P(z, R) 



-00 

J^(x) «ixp /tt Gx”^ erf (x/2G)] dx • 

•'o 


( 21 ) 


( 22 ) 


The power of the unscattered beam at an optical depth of T is, of course, 

—X 

F e . Thus, the presence of forward scattering has increased the detected 
■ o 

power by the factor 


A(t , G) 
s 


J (X) exp 

•'o 



A G X 


erf (X/2G) 


dX . 


( 23 ) 


A, which we call the amplification factor, is a function of two parameters; 

T , the scattering optical thickness, and G, the geometry factor, 
s 
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Finally, we may obtain the beam spread by substituting from Eqs . (11) and 

(16) into Eq. (14) : 

2 2 2 2 ^2 -2 

<r > = z /3a + z /3 + Y • (24) 

In general, the first term should dominate, except perhaps close to the 
point of entry into the medi\am. 

4.2 Non-Gaussian Phase Functions 

Although the Gaussian form in Eq. (17) is a popular model for the 
forward peak of the phase function, it is often a good idea to examine other 
models, to make sure that none of the results are simply an artifact of the 
Gaussian model. In this section, therefore, we shall examine a number of 
other functional forms which may be (and have also been) used to model 
anisotropic phase functions. We shall follow essentially the same steps as 
in the previous section, and present only the results, unless further 
explanation is necessary. 

4.2.1 Exponential Phase Functions 


i) 

P0|^) 

2 -a il; 

= a e V27T 

(25a) 


P(C) 

= a^(a^ + ^2^ -3/2 ^ 

(25b) 

where 

Q 

o 

= (1 + y^/6)"^^^ 

(25c) 


y 

= z n/a = X/ G • 

(26) 
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Thus 

and 

ii) 


Hence 


Note that, although Eq. (27a) implies' P(0) = the inclusion of the correct 
solid angle factor leads to a finite result for the amount of light scattered 
through any angle. In fact, Eq. (27a) has been employed by Bravo-Zhivotovskiy 
et al . (Ref. 6) to model the phase function of sea water. 

4.2.2 Binomial Phase Functions 

This time, we consider phase functions based on the functional form 

2,2 -y -1 

(1 + a ) . We will need the result that 

r°° 2 2 -11-1 11 2 

J (ni|i) (1 + a > ) ^ iji dip = (n/2a)^ K, (n/oc) / a r(y + l) (28) 

o y 

0 

where K is the modified Bessel function of the second kind. Thus if 

y 


A(T , G) = 
s 


.00 


0 


Jj^(X) exp 


T (1 + xVe dx 


<r^> = 2 T z^/a^ + zV3^ + Y ^ 

s 


P(i|j) = a 41 "^ e"“ ^/2T^ 


(25d) 

(25e) 

(27a) 


P(C) = a(a^ + / 2TT 


(27b) 


a 


= T y ^ In [ y//2 + (1 + |- y^) 


A(T , G) = 
s 


J^(X) exp|i^ /I Gy ^ In 


X/g/2 + il + ~x^G 


dx 


(27c) 

(27d) 


<r^> = 2 T / 3 / 3^ + Y~^ . 
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p^(4j) = 2 y (1 + ol\^) ^ ^ / 27t 

^ (?) = (?/2a)^ K (?/Q{) AT (y) 

y y 

a = /-rr r (p + i) [k, (y‘) l,, , (y") + K (y') L^Cy’)] / r(y) 
o s 2 *■ y y-j. y JL y 

where y' = yv^'-l 


and L is the modified 

y 

The expression for A 
obtain the beam spread: 


Struve function of order y. 
may be easily written down. 


For H > 1, we may 


<r^ = 


z^/3 a^(y - 1) 


2 2 

+ Z /3 + 


-2 


Note that if y is an odd half -integer , Eq. (29c) may be expressed in terms 
of exponential functions. For example, for y = 3/2 we find 


^ (3/2) = T [2/2y^-e (1 + 2 /2 y ^)1 

o s ^ 

A(T , G) = f” J, (X) exp { T [ 2/2 Gx ^ “ e ^'^^’^(1 + 2 A Gx } dx 

s 1 ® 


(29a) 


(29b) 


(29c) 


(29d) 


(30) 


( 30 ') 
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4.3 Beam Profile 


The spreading of a laser beam in a forward-scattering medixmi can be 

2 

most simply described via the beam spread parameter, <r >, which we have 
derived above. However, the profile of the expanding beam is also of interest 
and will now be considered. 

The formal expression for irradiance versus distance from the axis is 
given by Eq. (11'). In most cases, this expression is well-behaved. However, 
in the special case of 3r Y Eq* (11') will diverge. An alternative 

expression for N(r) may be obtained either by an integration by parts of Eq. (11') 
or by differentiating Eq. (12) . Adopting the second approach and setting 

/N -2 

I = F (2tt) , we obtain 

o o 


N(r) 


1 dP 
27TR dR 


1 2 -T I 

= - (2tt)“ F r“ e J, (X) e ° y dx (31) 

o Jq 1 o 

where the prime denotes differentiation of with respect to y (Eq. 26) . With 

the exception of the sea-water phase function, goes as y ^ for large y, and 
^ -2 

so goes as y , and convergence is assured. 


5 . 0- APPROXIMATE SOLUTIONS 

In an effort to simplify the above analysis, several authors have 
employed a number of approximations. In this section we shall examine two 
of these approximations, neither of which appears to be particularly useful 
in our problem. 
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5 . 1 Dolin-Fante Method 

Dolin (Ref. 8) and Fante (Ref. 9) have argued that the angular shape of 
the scattered intensity should be a much more slowly varying function than 
P(ip), and have thus extracted it from the integral on the right side of Eq. (3). 
After separating the scattered intensity from the unscattered, and Fourier 
transforming, they arrive at the following expressions 


^u “Gz 

I (z, ]]/ 0) = e Z n) 


(32a) 


and I (z, n/ O) = I (n, z n) 

~ ~ O ~ J 


dz ' exp [- G z ' - 


X(t, rif z n) dt] 


(jO G P [n (z - z ' ) ] 
o 


(32b) 


1 2 2 

where A(t, n, zn)=-r(iicriD“ |z r\ -tnl +(1-w)g 

4 O ' - - ' ^ 


(33a) 


A ( t ,ri , zH ) dt 




“ z')^ + (1 - 03 ) G (z - z') 
' o 


(33b) 


2 

and is defined by Eq. (17a) . 

Equations (32) and (33) may now be inserted in Eqs. (11) or (12) as required. 

>S 

Although P is no longer exponentiated, this result is complicated by the 
additional (finite) integration over z*. In the case of a Gaussian phase 
function, it is possible to reverse the orders of these two integrals, and 
perform that over n, to give 


A = e 


^ f r 2 1 : 

i ^ - 03 T dt exp jo3 T t - G / (— 03 T t 

O Jo o 30 


+ t^)] . 


(34) 
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Our calculations show that this approximation is reasonably accurate, 
except in those situations where T is large and G is small. In section 7 
(Numerical Results) , we will compare the predictions of Eq (34) with those 
of Eq. (23) . As it is not possible to perform any of the integrals for 
any of the other phase function models, we have limited our examination 
of this approximation to the case of the Gaussian phase function. 

5 . 2 Arnush-Stotts Method 

In order to extract analytic answers, Arnush (Ref. 10) and Stotts 

(Ref. 11, 12) have expanded P to second order before performing the integration 

to obtain (Series expansion of would yield the same result.) Arnush 

has used Bravo-Zhivotovskiy ' s (Ref. 6) sea water phase function, Eq. (27a), 

while Stotts originally used a Gaussian phase function, Eq. (17)^ and more 

recently the sea water phase function. This approximation is sufficient 

2 

to provide the correct values for both p(z,«») , and <r >. 

We start by re-writing the definition of as follows 

-1 

^ = 2tt 03 az (r)z) ^(t) dt 

° ° Jo 

-1 

= 2tt 03 az (nz) J (t\|3) P(l/3) 4; dl|3 dt . 

•'0 ■'o 

Expanding as a power series leads to 

2 2 ~2 

= 03 azd-qz 3i3 /12 + . . .) 
where xs defined by Eq. (17a) . 


(35) 


( 36 ) 



r 
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It is complicated, but reasonably straightforward to obtain the following 
expression for the beam spread parameter: 

<r^> = ^ 0 ) a + z^/3^ + (37) 

3 O 

Ignoring the higher order terms in Eq. C361 , we may insert this expression 
into Eq. (11^) to obtain 


N(z, r) = Fq ®xp (1 - 


T - r^/<r^>] / TT <r^> . 


(38) 


Similarly, integration of Eq. (12) leads to 


P(z, R) = F e 
o 


-d - W )T 
o 


1 - e 


2 2 ”I 

-R /<r >' 


(39) 


Ignoring 3 and y, we may re-express Eq. (39) in more familiar terms 

— (1 — Ui )T 1 — 2/ 

P(z, R) = F_ e ° ^ ^ 


(39') 


i.e., A(t^, G) 


. [. ■ 


exp (- 3 / T ) 

s 


(40) 


E:q>ansion of to second order in y is equivalent to an asymptotic 
o 

expansion to second order, in G or R Thus we may expect this approximation 

to be accurate for large values of G or R. However, its behaviour for small 
values of these parameters is quite different from that of the exact results 
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quoted in previous sections . Thus we cannot expect this approximation to 

prove particularly useful for our problem, as shown later by numerical comparisons. 

In fact, one finds values of A which are less than unity I 

For our problem, of course, we are concerned with small values of R and G, 

and hence we are interested in the behayiour of for large values of y, i.e.r 

o 

its asymptotic expansion. The phase function model of Eq. (27a) does not have 
an asymptotic expansion, due to the fact that P(0) = c». For the other cases 
we may easily show that 

Q ^ r q/y (41) 

o s 


where 


q = 2 tt 


.00 

P(i|;) dij;/a • 

0 


(42) 


Thus for a Gaussian (Eq. (17)), q = /tT; Eq. (25a) gives q = 1; and Eq. (29a) 

11 3 

gives q = y B (— , y + ”) , where B is the beta function. (For y = — , for example, 

q = 2 . ) 

For large values of y, the Arnush-Stotts approximation to goes to 
(minus) infinity, and so we cannot expect this approximation to accurately 
predict the power received by a small detector. 


6.0 EXACT METHOD OF TAM AND ZARDECKI 
The method of Tam and Zardecki (Ref. 7 ) is exact, at least in principle, 
but requires the evaluation of multidimensional intecrrals, the order of which 
is equal to the order of multiple scattering involved. We will restrict this 
discussion to the case of the Gaussian phase function only. 
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The Tam and Zardecki method consists in expanding exp in a Taylor 
series, before performing the integration over z' (Eq. • Thus, inserting 
Eq. (18) in Eq. (S'), and performing the Taylor expansion, yields 


m 


exp (^^ ) = 1 + E_ 

o m=l m 


m: 




dz. 


2 p 2 

dz^ exp { -n z^ / 4a } . (43) 


We may now perform the inverse Fourier transform (Eq. 11') to yield 


m 


-T 2 s 

N(z, r) = F e (tt z ) T N (z, r) 

o m=0 m I m 


(44) 


where N 


2-22 2-22 

= ^ s T /. ,_r B T, , } 

2 2 ^2 ^ ^22 ^ 2 ^ 
zy+B zy+B 


(45) 


and N 


m 


rl 

0 ^0 


[ dz ■ ■ * dz A ^ exp f-r^/z^ A 1 
J- 1 mm L / 


(45b) 


where A 


m 


- 2^2 --2 -2 -2 
a .L.z. +B +z y 
1=1 1 


(45c) 


We may note in particular that may be evaluated analytically in terms 
of the error function. The resulting expression is quite complicated, except 
in the case where B and y go to infinity, in which case we get 


Nj^(z, r) 


2 


a 


/tt [l - erf(g)] / 2g 


(46a) 


where 


g = r a/z . 


(46b) 
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Turning our attention to the power received, we obtain an e^^ression 
similar to Eqs. (44) and (45) , viz. 


.m 


— T S 

P(z, R) = F e Z — ~ P (z, R) 
o m=o m' I m 


(47) 


where P = 1 - exp 


2 

R -g . . T -} 


(48a) 


and 


m 


dz. 


o o 

dz {l - exp [“R /z A 1 } 

m m *' 


(48b) 


Note that, from Eqs. (47) and (48) , Eq. (13) may be obtained trivially. 

As with N^, P^ is also analytic , and in the simple case of B, Y we obtain 


^^(z, R) = 1 - e"^ + G /it [l - erf (G)] . 


(49) 


The number of terms required for the convergence of the series in Eqs. (44) 
and (47) grows steadily with T^, and so in some cases for large optical thicknesses 
it may become prohibitively expensive to use it. Nevertheless, these results 
have one use in that Tam and Zardecki (Ref. 13) have shown that the mth order 
terms in Eqs. (44) and (47) correspond to the contribution from mth order 
scattering. This in itself is a useful result. 

Another use suggests itself, however- The Gaussian phase function is simply 
a model, with the parameters ot and w available for adjustment to match "real" 
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scattering patterns. Since we now have a simple expression for the singly 
scattered contribution, we may compare it with that produced by a real phase 
ftmction and adjust cx and accordingly. Then we may use the results outlined 
above to estimate the multiply scattered contribution from such a phase 
function. Comparisons with second and higher order contributions are also 
possible. This should increase our confidence in the worth of results obtained 
from a model phase function. 


7 . 0 NUMERICAL RESULTS 

In this section, we shall present some typical results based on our exact 
formulation and the Arnush-Stotts type approximate method from selected computa- 
tional results. We shall examine 4 phase function models: Gaussian (Eq. 17), 

both exponential models (Eqs. 25a and 27a) , and the binomial model with \i = 3/2 
(Eq. 29a) . To simplify discussion, we shall refer to the phase function model 
of Eq. (25) as the exponential model, and that of Eq. (27) as the sea-water model. 

We start by examining the phase functions themselves. It is, of course, 
much simpler to plot the normalized phase function, P, rather than P, where P 
is defined by 

P = 2tt P / . (50) 

Unlike P, P is now a function of only one variable, aij>. In Figure 1, we plot 
P against a^, for 0 - aip - 3.5. 
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The next function we examine graphically is In Fig. 2 we plot 

/t against y for all 4 phase functions, as well as the Arnush-Stotts 
o s 

approximation. From this log-log plot, the asymptotic behavior of the 4 
phase functions is apparent, particularly that of the sea-water phase function, 
which has no asymptotic expansion. Although it is not obvious from this figure, 
the curve for the binomial phase function actually lies slightly above the 
sea-water curve for y values less than about 3. Finally, we note that for y 
values greater than 2, the results obtained by the Arnush-Stotts approximation 
differ markedly from those by our exact formulation, rapidly approaching large 
negative values for y greater than 4. 

We now turn to a discussion of the amplification factor. A, and the 
power received, P, as predicted by these 4 models, and also the Arnush-Stotts 
approximation. We have evaluated both A and P for G between 0.01 and 1.0, and 
between 0.5 and 15.0. (Throughout, we have assumed 3/ Y 

In the Appendix to this report, we have included a listing of the FORTRAN 
program used to generate this data, along with a brief explanation and sample 
output . 

In Figure 3 we plot A against G for a series of values, for the 
Gaussian model. In Figure 4 we plot P against (assuming unit incident 
power) for a series of values of G, again for the Gaussian model. Also shown 
on this plot is the transmission, T, which represents the power that would 
be received if all scattered light was lost. These two graphs clearly 
indicate the important role that forward scattering can play in the detection 
of transmitted beams, especially for optical thicknesses of the order of 


10 or higher. 
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In Figure 5 we plot A as a function of G for t =4.0 and 10.0, for 

s 

all 4 model phase functions, for our foirmulation as well as the Arnush-Stotts 

approximation. In the latter case, one finds values of A which are less than 

unity. Note that the binomial and sea water curves cross for both T values 

s 

(cf . Fig. 2) . For large values of G, we see that there is little to choose 
between the four phase function models. 

In Figure 6, we plot (A - 1.0) against G in log- log form for T = 1.0, in 
order to emphasize the linear relationship implied by Eq. (A3) in the Appendix. 

We see that for G less than 0.3, the integral term in Eq. (A3) makes a negligible 
contribution. In Figure 7, we plot (A - 1.0) against G for T = 5.0. Here we 
see that the integral term in Eq. (A3) is starting to make a contribution. Also 
in this graph we have included the Arnush-Stotts and Dolin-Fante approximation 
results. The Dolin-Fante result was not included in Figure 6 as it could not 
be distinguished from the exact result for the Gaussian phase function. 

8.0 CONCLUDING REMARKS 

The propagation of a laser beam in an optically dense medium such 
as a fog, dust stozm, or smoke is a problem of growing importance, both for 
communication and detection purposes. Although such dense media lead to a 
significant attenuation of the primary beam, much of the scattered radiation 
may still be found close to the beam axis and will, thus, be available for 
detection by a suitable detector , 

In this report, we have examined the spreading of a laser beam using the 
small-angle scattering approximation to the equation of transfer. This 
approximation appears eminently suited for the study of beam propagation in 
fog, dust, or smoke media, where the scattering phase function is highly 
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anisotropic. As well as the standard Gaussian model phase function, three 
other model phase functions have also been examined. The Gaussian functional 
form was used to describe the initial beam spread and profile, although the 
analysis is somewhat simpler (and the resulting expressions tidier) if the 
limiting case is taken. 

All the numerical results presented in this paper have been based on 
the assumption of a narrow, collimated beam. We may remark, however, that 
the results and expressions presented in this report (e.g., Eq. 20) may be 
applied with full generality. 

We have also examined a number of approximations which have been used to 

further simplify the expressions we have derived. The Arnush-Stotts 

approximation is quite suitable for use in the asymptotic regions, at large 

distances from the beam axis. However, the behavior of the solutions close to 

the beam axis is governed by the parameter q, the zeroth moment of the phase 
function. This moment weights the contribution from scattering through very 

2 

small angles far more highly than does the parameter ij; , the rms scattering 
angle, or third moment. In fact, for small R (i.e., small G) , one may 
expand Eq. (23) to first order (cf. Eq. A3) 

A(T^,G) ^ 1 + q G + . . . ( 51 ) 

Finally, we may remark that the method proposed by Tam and Zardecki makes 
a useful contribution by providing a connection between real (Mie) phase 
functions , and the parameters which must be used in the model phase functions 
used in this report. Further work on the applications of this method is 


recommended. 
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APPENDIX: COMPUTATIONAL DETAILS 

The numerical results presented in Figs. 3 to 7 were obtained from a 
relatively simple computer program, consisting of less than 100 executable 
statements. A full listing, and partial output, are included in this 
appendix. 

The task of this program is to evaluate the numerical integrals in 

Eqs.' (23) , (25d) , (27d) and (30') , for a series of values of T ^ between 

0.5 and 15.0, and a series of values of G up to 1.0. For comparison, 

the Arnush-Stotts approximation, Eq. (39) , is also computed. We have 

included the full results for x values of 4.0 and 10.0, which mav be 

s 

read in conjunction with Fig. 5. 

Although infinite integrals of this type are often handled by 
Gauss-Laguerre quadrature, this method was found wanting, due to the 
oscillatory nature of the integrand. Instead we have employed Simpson's 
rule, up to a finite cut off, allowing for the remainder of the integral 
by the following result: 

If 

f (X) - <P (X) for X ^ X 

and 

0(X) dx = $ is known, 

•*0 


then 



dX 




$ + 


•'0 


[f(X) - ^(X)] dx . 


(Al) 
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To apply this result, we note Eq. (41) , and choose X such that 

e ° 1 + q G/X , X 1 X . (A2) 


Thus 


(X) e dx ~ 1 + Tg q G + 


J (X) (e ° - 1 - q G/x) dx. (A3) 
± s 


Due to the wide range of values of G which we have used (three orders 
of magnitude) it is necessary to vary the step size accordingly. Thus we 
have used a step size equal to G, up to G = 0.22. A step size larger than 
this is unwise, due to the variation in the term. Thus, when we reach 
G = 0.22, a larger set of values is computed and stored, to be used for 
the remaining values of G. 

As pointed out above Eq, (41) , q is undefined for the sea water phase 
function, so we have set q = 0 in this case. As a result, we are forced 
to choose a considerably higher value of X in order to satisfy Eq. (A2) . 

In fact, if the sea water phase function was dropped from consideration, 
the time (and cost) of these calculations would be cut at least in half. 

As it is, the results we have obtained for this phase function must be 
considered distinctly less accurate than the others. 
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DIMENSION 83(2560)^ SIGNAL(4)#AMP(4 )»AS (3) >SA( 3)>Q( 3 I/GGI20) 
DATA SP>DX / 1.772453a5»0.05 / iQ / 12.0>2.0i6.0 / 
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X— OX/2.0 
C 

C SET UP BS ARRAY OF J1 BESSEL FUNCTIONS 

C 
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X«X+DX 
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5 CONTINUE 
C 
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C 
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C 
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C 
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DO 20 1-1/2560 
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IF(Y.LT.20) OMEGA-OMEGA-TAU+EXP(-Y)*( 1.0+2.0/Y) 

AMP ( AJ - AMP(4 )+Bl*EXP( OMEGA)*DX 
20 CONTINUE 
C 

C NOW CALCULATE THE APNUSH-STOTTS APPROXIMATION 

C (RESULTS FOR PHASE FUNCTIONS 3 C A ARE IDENTICAL) 

C 

00 25 1-1,3 

P-0 (I J*G+G/ A. 0/T AU 

ASd)-l.O-EXP(-P) 
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SIGNAL( I)-AHP(I)4TRANS 
25 CONTINUE 
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50 CONTINUE 
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1 ^SECOND LINE IS THE RECEIVED POWER (FOR UNIT INCIDENT POWER)*) 
70 CONTINUE 
STOP 
END 
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NORMALIZED PHASE FUNCTION, Zw P/a 



FIGURE 1. Normalized phase function P vs aij; for four model phase 
functions . 
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FIGURE 6. Subtracted amplification factor (A - 1.0) vs. geometry 


factor ( G) for T = 1 .0 . 
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